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Recent advances in computer technology have made it possible to use 
large Unite element models for shock and vibration analyses. One type 
of dynamic analysis is the calculation of responses and loads in the 
frequency domain for steady-state operating conditions. If substantial 
nonproportional damping levels are generated by energy dissipative 
components, such as bearings or hydraulic cylinders, then the calcu- 
lation of accurate frequency domain results can be computationally 
intensive for large models. To reduce the computational effort, it 
is common to assume that damping is proportioned to the mass and 
stiffness of the system. This proportional damping approximation can 
lead to significant errors in the frequency domain results. A method 
has been developed to produce very accurate results for this type of 
model without the large computational burden of a traditional non- 
proportional damping analysis. 


INTRODUCTION 

Recent advances in computer technology have made it possible, but not always practical, to use 
large finite element models for shock and vibration analyses. At one time, refined models were 
used only for linear static stress and deflection analyses which are less computationally demanding. 
Dynamic analyses were generally limited to smaller lumped parameter or coarse finite element 
models. However, today’s more aggressive design goals are promoting lighter weight structures 
that must operate effectively in higher performance environments. To help achieve these goals, 
dynamic analyses of refined finite element models are becoming a more accepted part of the design 
process. 

One type of dynamic analysis is the calculation of vibration responses and loads in the fre- 
quency domain for steady-state operating conditions. When damping levels are low and energy 
dissipation is well distributed throughout the system, proportional or modal damping approxima- 
tions usually produce sufficiently accurate results with a relatively modest computational effort. 
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On the other hand, if substantial levels of nonproportional damping are generated by energy 
dissipative components, such as bearings or hydraulic cylinders, then the calculation of accurate 
frequency domain results can be computationally intensive for large models. To reduce this effort, 
it is common to assume that damping is proportional to the mass and stiffness of the system. 
However, in this case, the proportional damping approximation can produce significant errors in 
the frequency domain results. A method has been developed to calculate very accurate results 
for this type of model without the large computational burden of a traditional nonproportional 
damping analysis. 


THEORETICAL BACKGROUND 

Finite element models can be used to represent the vibration behavior of structural systems in a 
steady-state operating condition. In matrix form, the system of N equations of motion can be 
written as follows: 


(-o> ! (M) + jw\C) + j\H] + [ir]) Wxff = {*■}„„, ( 1 ) 

where [M], [C\, [If], [K] are the physical mass, viscous damping, hysteretic damping, and stiffness 
matrices, respectively and {a} , {F} are the physical displacement and force vectors, respectively. 
The variable u> is the frequency (rad/ sec) and j = yj — 1 . 

When the finite element model represented by Eq. ( 1 ) exceeds several hundred degrees of 
freedom, the system is usually reduced prior to direct frequency domain calculations. A popular 
reduction method is to use M real normal modes of the system that span the frequency range of 
interest to create a transformation matrix and constraint relationship as follows: 


{ w };vxi — Mjvxm {^}mxi (^) 

where M <C N and the M columns of [it] are the real normal mode shape vectors and { 7 } is the 
modal displacement vector. By substituting Eq. ( 2 ) into Eq. ( 1 ) and premultiplying both sides 
of the resultant equation by [ , $ f ]' r , we have: 
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diagonal modal mass matrix 
modal viscous damping matrix 
modal hysteretic damping matrix 
diagonal modal stiffness matrix 
modal force vector 


(3) 


The reduced system represented by Eq. (3) is an approximation to the original system in Eq. 
(1). The modal displacement vector { 7 } becomes the new set of independent coordinates and 
the original physical displacement vector {u} is back-calculated through the mode shape matrix 
in Eq. ( 2 ). The modal displacements in Eq. (3) can be calculated at each frequency of interest 
through a frequency-dependent matrix inversion in this way: 
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{ 7 } = (-u.’l'm.J + J»[c]+,-[k] + r'fcj)'' {/} (4) 

In the reduced system of Eq. (3) and Eq. (4), the modal mass and stiffness matrices are diagonal, 
but the damping matrices can have large off-diagonal terms requiring a fully populated complex 
matrix inversion for each solution frequency. Although the system has been reduced, a large 
number of retained modes and/or large number of frequency steps can produce a computationally 
intensive solution. If a proportional damping approximation is made at this stage, the off-diagonal 
terms in the damping matrices are ignored and the resultant equations of motion are then fully 
uncoupled. In this case, the frequency domain solution in Eq. (4) becomes very efficient because 
a scalar, rather than matrix, inversion is required at each frequency of interest. Neglecting 
the effects of these off-diagonal coupling terms in the damping matrices, however, can generate 
substantial errors in the frequency response and load calculations. 

The presence of discrete damper components in the model does not usually influence all of the 
system modes to a significant extent. Only those modes that have a substantial amount of relative 
motion across the dampers will be strongly affected by their energy dissipative properties. When 
little relative motion exists across the dampers for a given mode, the dampers are not effectively 
exercised and little energy dissipation is produced for that mode. Therefore, Eq. (3) can be 
partitioned so that modal coordinates corresponding to the modes that significantly exercise 
damper components are separated from the other modal coordinates in this way: 



where the subscripts n, p denote the ‘nonproportionally’ and ‘proportionally’ damped partitions, 
respectively. 

The new analysis technique presented in this paper is a hybrid of the traditional propor- 
tional and nonproportional damping solution methods. The off-diagonal terms in the c pp and h pp 
quadrants of the damping matrices in Eq. (5) should be small compared to the diagonal terms of 
their respective quadrants. Ignoring these off-diagonal terms and replacing those quadrants with 
diagonal or proportionally damped matrices is a small approximation. The larger off-diagonal 
terms in the c nn , c np , c pn , /i nn , h np , h pn quadrants are retained to represent the significant non- 
proportional damping effects. This approximation becomes the basis for a condensation of {7 P } 
onto the {7 n } modal coordinates which eventually creates a more efficient method of solving Eq. 
(3) with a small loss of accuracy. The concept of partitioning and condensing is similar to the 
Guyan reduction technique [1] although the present application is very different. 

The new hybrid formulation begins with the expansion of Eq. (5) for the pn and pp quadrants 
as follows: 

(-u> 2 r'mppj + jo: ["C PP J + j [ "hpp^l + [^ PP ^j) {7p} = {/p} - (jw [ c p«] + 3 [ h pn ]) { 7 -} (6) 

where [~~'c PP -J , [~"h PP -J are the diagonal (proportionally damped) approximations to the original 
damping matrix quadrants. From Eq. (6), the vector {7 P } can be expressed in terms as {7n} as 
follows: 
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where 


{.Ip) — [ A>P--] ({/p} U<*> [Cpn] + J [^pn]){7n}) 


( 7 ) 


[ Dpp ^ j — ( ^ + j<*> r-cp ^ + j r-vj [ ^w-j ) 

By substituting Eq. (7) into the expansion of the nn and np quadrants of Eq. (5), we have: 

[Am] {7n} + (jw [Cnp] + j [J*np]) ["ApJ” 1 ({ f p } - (joj [c pn ] + j [hpn]) {7n}) = {/n} (8) 

where 

[Am] = (-u> 2 [ m nn ] + ;u> [c nn ] + j [fi„„] + [ k nn ]) 

Eq. (8) can be simplified and inverted to solve for { 7 ,,} as follows: 

{7n}= [P»n]"'{/n} (9) 

where 

[Am] = ([Am] + U 1 [C„ p ] rAp ^" 1 [ C Pn] 

+ w [ c np] [ -Dpp~J [hpn] + tV [hfip] [ Dpp — ] [Cpn] 

+ iM r^vr 1 i u) 

and 

{/.} = {/»> - w + ; ia»pD r^v-T' {/,} 

The calculation of {7 n } through the hybrid solution in Eq. (9) provides an accurate, but 
more efficient, alternative to the traditional nonproportional damping solution in Eq. (4). The 
hybrid system matrix, [^Xm] although fully populated, is usually significantly smaller than the 
modal system matrix in Eq. (4). Because the number of calculations required for the inversion 
of fully populated matrices increases cubically with the matrix order, substantial savings can be 
gained through order reductions. Additional effort, of course, is required for the calculation of 
|jD nn ] and j/ n j at each solution frequency and offsets some of the computational savings gained 
from the matrix inversion of a smaller system. However, the small approximation of replacing a 
nearly-diagonal [D pp ] with a diagonal p* D pp matrix quadrant makes these additional calculations 
for the generation of i) nn j quite manageable. 

Once { 7 „} has been determined for each frequency using Eq. (9), the remainder of the 
modal displacements { 7 P } can be calculated using Eq. (7). Finally, the physical displacements 
{u} of the original dynamic system can be calculated using Eq. (2) thus completing the frequency 
domain solution of Eq. (1). 

APPLICATION 

Recently, the authors performed an evaluation of a rotating equipment design intended for 
marine application. This evaluation included a frequency response analysis using a very large 
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finite element model of the system where more than 200 resonant modes were needed to span 
the frequency range of interest. The solution was to be calculated at more than 1500 spectral 
lines to provide the desired resolution over the frequency range of interest. Because significant 
levels of nonproportional damping were present in the model and accurate results were required, 
a traditional nonproportional damping analysis was attempted. Using MSC/NASTRAN software 
on a CRAY mainframe computer system, a partial solution over a narrow frequency range was 
performed. The CPU time used for each frequency step of solution was linearly extrapolated to 
estimate the time for a complete solution. The required CRAY CPU time was estimated to exceed 
40,000 seconds and the solution was judged to be impractical using conventional approaches. As 
a result, the alternative method presented in this paper was developed. 

The new hybrid method was first applied to another large finite element model whose size 
and construction was similar to the model previously discussed. The 71,000 degree of freedom 
model was divided into six substructures where each substructure was generated using the Craig- 
Bampton formulation [2]. The residual system of assembled substructures consisted of 4,300 
degrees of freedom. It was dynamically reduced prior to frequency response solution using 123 
real normal modes that spanned the frequency range of interest as outlined in Eqs. (1), (2), and 
(3). Viscous damper and hysteretic damper elements were used to represent the energy dissipative 
effects of fluid-film bearings and elastomeric mounts, respectively. The modal damping matrices 
were nearly fully populated with large off-diagonal terms indicating that significant levels of 
nonproportional damping were present in the system model. 

The frequency response solution of the reduced system with 123 modal coordinates was 
calculated at 400 spectral lines for each of the three methods discussed in this paper which 
include: 

• Method 1: Traditional nonproportional damping solution 

• Method 2: Traditional proportional damping solution 

• Method 3: New hybrid solution 

Figure 1 presents a comparison of calculated frequency responses for each of the three 
solution methods evaluated at the same structural location. Method 1 results are ‘exact’ and 
differ substantially from the approximate results of Method 2. The computational time required 
for the accurate results of Method 1, however, was 52 times greater than the approximate Method 
2 solution. On the other hand, Method 3 results from the hybrid solution compare favorably with 
the ‘exact’ solution from Method 1, but required only 8% of the computational time needed for 
that solution. 

For this particular hybrid solution, the reduced system matrix in Eq. (5) was partitioned 
into 31 ‘nonproportionally-damped’ and 92 ‘proportionally-damped’ modal degrees of freedom. 
In other words, 31 real modes were judged to be significantly influenced by the discrete damper 
elements in the model. The remaining 92 real modes were not significantly affected and were 
subsequently approximated with 92 ‘proportionally-damped’ modes, thus producing an efficient 
solution with little sacrifice in accuracy. 

The matrix partitioning process outlined in Eq. (5) was not a straight forward task. As 
more modes are included in the ‘nonproportionally-damped’ partition, the accuracy as well as the 
computational expense will increase and eventually converge to the traditional nonproportional 
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damping solution. Several partitioning schemes were developed and studied based on the relative 
magnitudes of off-diagonal to diagonal terms in the damping matrices. In the future, more 
sophisticated partitioning algorithms are expected to generate even more economical solutions. 


ACCELERATION (dB> 



Method 2 : Proportional Damping Solution (CPU = IX) 
ooooooo Method 3: Hybrid Solution (CPU = 4X) 


Figure 1 

Comparison of Frequency Response Calculations for Different Solution Methods 


CONCLUSIONS 

A method has been presented to accurately and economically calculate steady-state frequency 
responses based on the analysis of large finite element models with nonproportional damping 
effects. The new method is a hybrid of the traditional nonproportional and proportional damping 
solution methods. It captures the advantages of each conventional approach without the burden 
of their respective shortcomings, as demonstrated with comparative analyses performed on a large 
finite element model. 
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